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Summary 

A thermal model which simulates combined conduc- 
tion and phase change characteristics of thermal energy 
storage (TES) material is presented. Both the model 
and results are provided herein for the purpose of 
benchmarking the conduction and phase change capa- 
bilities of recently developed and unvalidated microgravity 
TES computer programs. 

The model was originally developed to support the 
design and analysis of the NASA Lewis Research Cen- 
ter TES experiment. This report also serves as the 
preliminary documentation of the thermal models con- 
structed for this experimental effort. The phase change 
approach developed for the model can be applied to 
most commercial thermal analyzers and also can be 
used for other phase change applications. It is again 
emphasized that the principal intent of this report is to 
present a model and results for evaluating (benchmarking) 
computer codes that have been recently developed to 
predict the behavior of phase change materials in a 
microgravity environment. 

Specifically, operation of the TES-1 experiment is 
simulated. The heart of the TES-1 experiment is the 
TES specimen. The specimen is mounted in an approved 
shuttle payload container known as a Get-Away-Special 
(GAS), which is also included in the model. The TES 
specimen consists of the torus shaped TES canister 
which contains the phase change material (PCM), the 
conductor rod, and the radiator flare. A heater radiates 
heat to the outer radius of the TES canister during 
melting of the PCM. After the lithium fluoride PCM is 
completely melted, the heater is turned off and latent 
heat is transferred from the PCM to the conductor rod 
at the TES canister inner radius. The removed heat is 
radiated to the environment from the TES radiator flare 
during solidification. A series of five melting/freezing 
cycles are simulated. 

A two-dimensional SINDA85 model of the TES 
experiment in cylindrical coordinates was constructed. 
The phase change model accounts for latent heat stored 
in, or released from, a node undergoing melting and 
freezing. However, the volume change of the LiF dur- 
ing the phase change process is neglected. Other phe- 
nomena not considered in the model include buoyancy 
driven and Marangoni flows and the semitransparent 
behavior of the lithium fluoride (LiF). 


Three types of results provided for comparison with 
other models are transient temperature profiles, color 
temperature contour plots, and tabulated results. Tabu- 
lated results are peak temperatures, freezing times, melting 
times, and heat losses following each freeze and thaw 
period. Results from four simulations are presented. 
Three simulations were performed to investigate the 
effects of the GAS lid mass on the melting and freezing 
times of the PCM. The fourth simulation uses a simple 
approximation to account for combined radiation and 
conduction transport through the liquid LiF (the 
semitransparent nature of the solid phase was not included). 
Results from the lid mass comparison indicate negli- 
gible effect on the freezing and melting times for GAS 
lid masses from 9.1 to 34.1 kg. The semitransparent 
approximation resulted in decreases of 3.9 percent in 
freezing times, 3.8 percent in melting times, and 20 K 
in peak canister temperatures with respect to conduc- 
tion only. 

Introduction 

Thermal energy storage (TES) is an attractive tech- 
nology for solar dynamic power systems proposed for 
space applications. A TES system would deliver a 
near-constant temperature working fluid to the turbine 
during a sun-shade orbit. A typical solar dynamic power 
cycle which includes the TES material, lithium fluoride 
(LiF), is shown in figure 1. Energy from the concen- 
trator is stored as latent heat in the LiF. One proposed 
receiver design utilizing thermal energy storage is 
illustrated in figure 2. The TES material or phase change 
material (PCM), is actually contained in small (1.27- to 
2.54-cm long) torus shaped canisters that surround the 
heat pipe. A principal feature of this configuration is 
that heat (solar flux) is input to the TES material from 
the outer TES canister radius, while heat is extracted 
by the heat pipe at the inner canister radius. 

Various fluoride salts (e.g., LiF, LiF-CaF 2 , NaF) which 
undergo phase change in the temperature range required 
for the proposed power cycle have been identified as 
candidate thermal energy storage materials (ref. 1). 
However, fluoride salts exhibit as much as a 30 percent 
change in volume during the phase change process. 
This leads to the formation of voids during the solidi- 
fication process. Two concerns which emerge as a 
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Turbine Compressor 

Figure 1 . — Typical solar dynamic power cycle. 



result of void formation are the potential for “ratcheting” 
of the TES containment wall and development of “hot 
spots” on the TES canister (containment). Ratcheting 
is the distortion that occurs if the liquid PCM is unable 
to expand into the void and instead creates stresses on 
the canister wall as the PCM continues to liquify dur- 
ing the melting process. Hot spots occur on the canis- 
ter wall in regions where voids exist between the wall 
and salt. The present lack of understanding of the 
freeze-thaw process in a microgravity environment 
coupled with large volume changes results in conser- 
vatively designed TES/receiver systems (refs. 1 and 2). 
The prediction of void behavior is also of interest to 
other space related efforts such as possible SP-100 cold 
restarts. 

Computer capabilities to predict void nucleation, move- 
ment, and final location/shape in a microgravity 
environment are being developed (refs. 2 to 4). The 
computer programs or models must obviously be evalu- 
ated to establish their credibility. Confidence in the 
model is initially established by verification with known 


analytical solutions. Additional capabilities can be assessed 
by benchmarking the unvalidated computer program 
with a validated computer program. The validated model 
does not necessarily provide an evaluation (benchmark) 
of all capabilities of the unvalidated model; thus, it can 
be less sophisticated than the unvalidated model. The 
final step in the evaluation process is the validation of 
the code with actual experimental data. 

Microgravity data to validate the microgravity TES 
models are currently nonexistent. However, a set of 
space flight experiments, the NASA Lewis TES flight 
experiments, are presently being designed and constructed 
to provide both insight into the void behavior and the 
data to evaluate these computer models (ref. 1). In 
support of the TES experimental design effort, ther- 
mal analysis has been performed that can also be used 
in the evaluation (benchmarking) of the microgravity 
TES computer programs. The model and results from 
this analysis are presented for the purpose of benchmarking 
some of the capabilities of the unvalidated TES com- 
puter models. 

Combined conduction and phase change within the 
phase change material (PCM) are simulated in the model 
presented. The PCM volume change during phase change 
is not simulated. Since the model was originally con- 
structed for design and analysis purposes, the assump- 
tions and modeling approach were strongly influenced 
by design considerations. However, the emphasis of 
this report is on presenting sufficient information for 
comparison with results from other models or computer 
programs. 

Symbols 

C thermal capacitance J/K 

Cp specific heat, J/kg-K 

Gy linear conductance between the i and j 

nodes 

G rad radiation conductance between the i and 

j nodes 

h contact conductance, W/m 2 -K 

h lat latent heat, J/kg 

K thermal conductivity, W/m-K 

k counter - kth iteration 

m i mass of ith node, kg 

N number of nodes 
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QERR convergence limit for phase change initiation 
or completion logic 

Q exs i excess latent heat of ith node at phase change 
completion time step 

Q i lat cumulative latent heat of the ith node 

Tj n temperature of the ith node at time n 

T- n temperature of the ith node at time n 

T lat solid/liquid phase change temperature, K 

t n discrete time value immediately prior to phase 

change completion 

t n+ j discrete time value immediately following phase 
change completion 

X i liquid mass fraction 

At time step 

m At m +i time step in which phase change initiation 

occurs 


results to other numerical or predicted results. Finally, 
a detailed discussion of the phase change approach is 
provided because the logic can easily be incorporated 
into other thermal models with different phase change 
applications. 

Model Geometry and Thermal Properties 

The thermal model geometry of the proposed TES 
design was constructed with the use of PATRAN. The 
PATRAN model was then converted to a SINDA model 
in cylindrical coordinates. Since PATRAN is a finite 
element computer program while SINDA is a lumped 
parameter (lumped volume) computer program, PATRAN 
nodes and SINDA nodes are not equivalent. The term 
“node” in this report shall always imply a SINDA node 
unless otherwise stated. The node diagrams for the 
overall model and TES specimen are shown in figures 3 
and 5, respectively, while dimensions for the overall 
model and TES specimen are found in figures 4 and 6, 
respectively. The model is two-dimensional in the radial 
and axial directions and is comprised of 127 nodes. 


n At n+1 time step in which phase change completion 
occurs 
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e emissivity 

e ef f effective multilayer insulation (MLI) emissivity 

T| latent heat content error 

Subscripts: 
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Thermal Model Description 

Reasonably complete descriptions of the model geometry, 
model properties, phase change logic, orbit simulation, 
and general modeling approach are provided in this 
section. A familiarity with the modeling approach and 
assumptions are important for the correct interpretation 
of the model results. In addition, an understanding of 
all the above items is critical for comparison of model 
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Figure 3. — GAS/TES node diagram. 
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Figure 4. — GAS/TES model dimensions (cm). 


Node distribution is as follows: 


TES specimen nodes 
LiF nodes 

108 

20 

MLI 

5 

Lid, GAS, and support structure 

12 

The node distribution reflects the level of detail that 
was used in representing various components. Hence, 
from the above node distribution, we observe a rela- 

tively detailed TES specimen model, 

while the GAS 


walls, GAS lid, and MLI are coarsely represented. When 
attempting to compare the results with other numerical 
analyses, it is also helpful to give a breakdown of 
diffusion nodes, arithmetic nodes, and boundary nodes. 
Diffusion nodes are nodes with finite capacitance val- 
ues (i.e., time varying nodes), while arithmetic nodes 
have zero capacitances (also referred to as massless 
nodes). Rigorously, the temperature of an arithmetic 
node at each time step is simply the thermal conduc- 
tance weighted average temperature of all nodes in 
contact with the given node. Boundary nodes are nodes 
with constant temperatures. The distribution of diffu- 
sion, arithmetic, and boundary nodes in the thermal 
model is, 70, 55, and 2, respectively. Most, but not all, 



Figure 5 — TES specimen node diagram. Circles denote arithmetic 
nodes. 



Figure 6.— TES specimen dimensions (cm). 
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of the linear conductances were generated from the 
PATRAN model, while radiation conductances were 
computed with the use of TRASYS and then were 
incorporated into the SINDA thermal model. The num- 
bers of linear, nonlinear (radiation), and total conduc- 
tances are 932, 202, and 1134, respectively. 

The TES specimen is mounted in an approved shuttle 
payload container previously known as a Get-Away- 
Special (GAS) and more recently designated as a Complex 
Autonomous Payload (CAP). The GAS (or CAP) is a 
cylindrical aluminum container with an internal diam- 
eter and length (user length) of 0.5017 and 0.7176 m, 
respectively. In figure7, the TES specimen is shown 
mounted in the upper portion of the GAS. The battery 
and electronics will be located in the region below the 
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Figure 7.— TES/GAS configuration. 

TES specimen. The surrounding MLI and mounting 
base assembly thermally isolates the TES specimen from 
the battery and electronics compartment and the GAS 
walls. The mounting base assembly consists of two 
silica spacers, four Haynes- 188 bolts, and the molybde- 
num heater bracket, which is sandwiched between the 
spacers. As illustrated in the GAS/TES node diagram 
(fig. 3), the Haynes-188 mounting bolts and silica spacer 


are combined. The thermal capacitances of the four 
bolts and silica spacer are summed together to obtain 
capacitances of the mounting base nodes, nodes 109 
and 9109. Thermal conductances of the mounting base 
nodes were computed by assuming parallel conduc- 
tances (resistances) between the bolts and the spacer. 
The axial thermal resistance of the thin molybdenum 
heater bracket is negligible relative to the spacer ther- 
mal resistance. However, on each side of the heater 
bracket a contact conductance, h, of 400 W/m-K was 
included in the thermal model. 



An enlarged view of the TES specimen is shown in 
figure 8. The conductor rod, radiator flare, and the 
TES canister are made of Haynes-188 alloy. For this 
configuration, the heat is radiated from the heater to 
the TES canister. Heat is transferred from the heater in 
the following ways: radiated to the canister wall, radi- 
ated through the upper and lower gaps between the 
canister wall and heater, and radiated from the heater to 
the MLI. Conduction from the heater along the heater 
mounting bracket is not included in the model. How- 
ever, the model is flexible enough to allow a bracket to 
be included in subsequent analyses. For assembly pur- 
poses, the radiator flare is detachable from the conduc- 
tor rod. Therefore, a radiator flare/conductor rod contact 
conductance of 400 W/m-K is also included in the model. 
While the MLI that surrounds the TES specimen is 
fixed, the MLI blanket, or cover, above the radiator 
flare is movable. A shutter mechanism is used to open 
the blanket and expose the radiator flare to the GAS lid 
during the solidification (cooling) portion of the cycle. 
During the melting (heating) portion of a cycle, the 
MLI cover is closed to insulate the radiator surface 
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from the surrounding environment. The simulation of 
the alternating heating and cooling portions of the orbital 
cycles is discussed later. A total heater power of 260 watts, 
or a heat flux of 15.1 kW/m 2 , was applied to the heater 
during the melting portion of each cycle. In each simu- 
lation, power was applied until all the PCM was melted. 

The thermophysical properties used in the model are 
summarized in tables I to IV. Temperature dependent 
specific heats and thermal conductivities associated with 
the model materials and components are given in table I. 


Constant values of densities and latent heat of fusion 
are shown in table II. Emissivities of various surfaces 
that were used to calculate radiation conductor values 
are given in table III. Contact conductances are given 
in table IV. Contact conductive values were chosen 
from a range of values provided in reference 6. 

The spacer thermal conductance values shown in table I 
refer to the mounting base spacer below the TES speci- 
men, as indicated in figure 7. The spacer thermal con- 
ductances were computed by assuming the thermal 


TABLE I.— TEMPERATURE DEPENDENT PROPERTIES 


[Specific heat, Cp, thermal conductivity, K, thermal conductance, U.] 


Specific heat for Haynes 188 alloy 2 








Temperature, K 

273 

673 

873 

1073 

1273 

1473 


Specific heat, J/kg-K 

398 

486 

523 

557 

590 

624 


Thermal conductivity for Haynes 188 alloy 3 








Temperature, K 

311 

644 

811 

1033 

1200 

1477 


Thermal conductivity, W/M-K 

10.8 

17.0 

19.9 

24.1 

25.8 

30.5 


Specific heat for LiF b 








Temperature, K 

500 

1121.3 

1121.4 

1500 




Specific heat, J/kg-K 

1631 

1631 

2453 

2453 




Thermal conductivity for LiF b 








Temperature, K 

500 

600 

800 

1000 

1121.2 

1121.3 

1500 

Specific heat, J/kg-K 

8.08 

6.61 

5.61 

5.98 

6.34 

1.73 

1.73 

Specific heat for aluminum (from ref. 5) 








Temperature, K 

200 

250 

300 

400 

500 



Specific heat, J/kg-K 

797 

859 

902 

949 

997 



Thermal conductivity for aluminum (from ref. 5) 








Temperature, K 

200 

250 

300 

400 

500 



Thermal conductivity, W/M-K 

237 

235 

237 

240 

236 



Thermal conductance for spacer (Node 109) 








Temperature, K 

311 

811 

1033 

1477 




Thermal conductance, U, W/K 

0.01376 

0.02133 

0.02482 

0.03014 





Property values obtained from Haynes Alloy No. 188 Manufacturer Physical Property Brochure, Haynes International Inc., Kokomo, 
IN. 

b Property values obtained from Hitchhiker Shuttle Payload of Opportunity Carrier Customer Accommodations and Requirements 
Specifications. HHG-730- 1503-05 (draft copy 2/90), NASA Goddard Space Flight Center, Greenbelt, MD. 


TABLE II.— CONSTANT PROPERTIES 


Density dcg/m 3 

Haynes 188 8980 

LiF 1690.3 

A1 at 300K (from ref. 5) 2701 

LiF latent heat, kJ/kg 1 037 


TABLE III— MODEL EMISSIVITIES 


Component 

Emissivity 

GAS lid outer surface 

0.75 

Gas lid inner surface 

.75 

GAS wall inner surface 

.75 

GAS wall outer surface 

.75 

MLI inner surface 

.10 

MLI effective emissivity 

.0025 

PCM void 

.5 

TES canister wall 

.7 

Heater inner surface 

.55 

Heater outer backing 

.5 

Radiator flare top 

.5 

Radiator underside 

.5 

Conductor rod neck 

.5 
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TABLE IV. -MODEL CONTACT 
CONDUCTANCES 


Interface 

Conductance, 


h, 


W/m 2 -K 

Radiator flare to conductor rod 

400 

Spacer to heater bracket 

1200 


resistance of the four Haynes- 188 bolts and the ceramic 
spacer to be in parallel. The thermal capacitance of the 
ceramic spacer relative to that of the bolts was neglected. 
Therefore, the total spacer thermal capacitance is the 
sum of the four bolt thermal capacitances. 

As depicted in figure 7, the MLI blanket that sur- 
rounds the TES specimen is represented by a series of 
arithmetic nodes. The inner surface, facing the TES 
specimen, is given an emissivity (£ = 0.1) typical of an 
MLI foil surface. An “effective” emissivity, £ eff , rep- 
resents the radiative resistance of the MLI nodes. The 
effective emissivity value is chosen to yield a desired 
or acceptable average heat loss through the MLI during 
a complete melting of the PCM. For example, the e eff 
of 0.0025 given in table III yields an acceptable MLI 
average heat loss of -52 W during PCM melting. 
Determination of the appropriate £ eff is an iterative 
procedure of guessing £ eff , running the melt/freeze 
simulation, and adjusting the new guess of £ eff based 
on prior MLI heat loss prediction. 

Most of the heat rejected during the TES experiment 
is radiated to the GAS lid. The lid is considered to be 
in excellent thermal contact with the rest of the GAS 
canister. Therefore, to minimize the temperature rise 
in both the lid and canister, it is necessary to reject as 
much heat as possible to the environment. According 
to a draft GAS requirements document, 1 various options 
are available for heat rejection from the GAS canisters. 
One such option calls for painting the outside surfaces 
of the GAS canister white to obtain an emissivity of 
0.86. Another option may be to paint the inside sur- 
faces of the GAS canister to obtain emissivities of 
approximately 0.8. As indicated in table III, a conser- 
vative emissivity value of 0.75 was applied to all GAS 
surfaces in the present thermal model. If GAS emis- 
sivities greater than 0.75 are obtained, then the heat 
rejection capability will improve. 

Phase Change Procedure and Logic 

The phase change modeling approach is described 
here in detail because this approach is different from 


hitchhiker Shuttle Payload of Opportunity Carrier 
Customer Accomodations and Requirements Specifica- 
tions. HHG-730-1503-05. (Draft copy 2/90), NASA 
Goddard Space Flight Center, Greenbelt, MD, pp. 2-44. 


the standard approach, which uses the effective spe- 
cific heat Cp eff (ref. 7), that can be applied to most 
thermal analyzers. One objective of the model is to 
provide results that will be used to benchmark com- 
puter codes currently being developed to predict the 
behavior of PCM (refs. 2 to 4). Therefore, a thorough 
understanding of the actual capabilities and limitations 
of the phase change logic is necessary to determine 
which characteristics of the unvalidated computer pro- 
gram are actually being benchmarked. The following 
discussion details the latent heat simulation and the 
initiation and completion of the phase change logic. 
The development of the latent heat simulation (phase 
change) logic is described in appendix A. 

The present phase change model accounts for latent 
heat stored in, or released from, a node undergoing 
melting and freezing. The phase change model cor- 
rectly accounts for latent heat as part of the nodal 
energy balance. However, other phenomena associated 
with the phase change of PCM’s are not considered. 
The model does not account for the volume change of 
the LiF during the phase change process. Buoyancy 
driven and Marangoni flows during the liquid phase are 
also not considered. Specifically, general conduction 
and solid-liquid phase change predictions from other 
computer codes can be benchmarked. 

Most commercial thermal analyzers do not explicitly 
include a latent heat term as part of the energy equation 
algorithm. However, as described below, it is relatively 
easy to add the capability of simulating nodal latent 
heat change. The most common approach is to increase 
the thermal capacitance of a node so that the sensible 
heat, Cp eff AT lat , over the temperature interval defined 
by AT| at is equal to h lat (ref. 7). Thermal analyzers 
that provide access to the basic solution variables (i.e., 
nodal temperature, thermal capacitance, and thermal 
conductance) allow a slightly different approach which 
is numerically simpler and more robust. This more 
robust approach has been incorporated into the current 
model. The major difference between the two approaches 
is that the infamous “machine roundoff error,” the chief 
disadvantage of the more common approach that uses 
effective specific heat Cp eff (ref. 7), is used advanta- 
geously to maintain a constant nodal phase change tem- 
perature in the current approach. In addition, a constant 
phase change temperature is not possible in the effec- 
tive Cp approach; rather, a small phase change AT is 
necessary (ref. 8). 

When phase change is initiated, the node thermal 
capacitance, C i? is set equal to a very large value (e.g., 
10 25 ). A sufficiently large value results in an 
extremely small temperature change that is rounded off 
by the machine to zero. Therefore, the node tempera- 
ture remains constant. Simulation of phase change is 
accomplished by cumulatively summing the net heat 
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flow into the node at each time step as given by 
equation (1). 



During a melting process and a solidification process, 
the liquid mass fraction of the ith phase change node 
is computed from equations (2) and (3), respectively. 


X; = ^ lal (melting) (2) 


Xj = 1 - Slj^L (freezing) (3) 

m ihlat 

When sufficient heat is accumulated or lost to com- 
pletely melt or solidify the node, then the liquid mass 
fraction, X if is equal to 1 or 0, respectively. The phase 
change process is then terminated by computing the 
appropriate value of C j required for sensible heat 
process. 

The explicit steps used to simulate melting and 
solidification for the present analysis are presented in 
the flowchart illustrations, figures 9 and 10, respec- 
tively. The primary function of this logic is to maintain 



Figure 9.— Melting logic. 
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Figure 10. — Solidification logic. 


the temperature constant and compute the liquid mass 
fraction during phase change. In addition, the state of 
the phase change nodes and the direction of phase change 
are tracked through the use of the flag, NPHASE. This 
is actually an array in the model that contains an ele- 
ment for each phase change node. For a given node, 
NPHASE values of 1, -1, 2, -2, 3, and 3, -3 indicate that 
the node is melting, solidifying, initiating completion 
of melting, initiating completion of solidification, liq- 
uid, and solid, respectively. The flowchart boxes in 
figures 9 and 10 that identify the initiation or comple- 
tion of phase change are expanded in figures 11 to 14. 


The question of when to apply the phase change 
logic must now be addressed. Since each time step, At, 
is finite, the probability is low that the temperature T i 
will be exactly equal to T lat at a given discrete time, 
t m . The discrete times that bound the unknown time at 
which melting of the ith node begins, t lat , are denoted 
as t m and t m+1 , while the time step between t m and 
t m +i is defined as m At m+1 . For melting, the times, t m 
and t m+1 can be determined by satisfying the criteria 
T i,m < T lat anc * T im+l > T laf For solidification, the 
above inequality signs are reversed. Therefore, the 
fundamental criterion that is used to identify when phase 
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Figure 11. — Melting initiation logic. 


change should begin is when the temperature differ- 
ence (T im - T lat ) changes sign from (Tj m+1 - T lat ). If 
the temperature difference (overshoot) between T m+1 
and T| at is small, then the sensible heat overshoot, 
Qi,os> is negligible, and t m+1 =t lat and T i m+1 = T lal . 

Qi,os = (T i>m+ i - T te ) (4) 

m&m+l 

The simplest approach to correct for non-negligible 
values of Qj os is to repeat the simulation using a 
smaller time step. When a limited number of simula- 
tions are required, this may prove to be the quickest 
approach. Disadvantages of this method are the increase 
in simulation time for the reduced time step and the 
additional number of simulations. An alternative is to 
develop a method that computes the initial latent heat 
change of the ith phase change node during the m At m+1 
time step. If the finite difference solution is explicit in 


time, then the latent heat contribution, Qj lat , during the 
m At m +i time step is simply equal to Qj os . As described 
in appendix A, no iteration is required for an explicit 
solution routine. However, in fully implicit or implicit/ 
explicit (e.g., Crank-Nicholson) finite difference solu- 
tions, all node temperatures are solved for simulta- 
neously at each given time, say t m+1 . Therefore, the 
latent heat contribution during the m At m+1 time step 
must be determined iteratively for implicit or implicit/ 
explicit solution routines. 

An iterative scheme was devised that corrects for the 
temperature overshoot by computing the latent heat 
contribution during the time step that phase change was 
initiated, m At m+1 . For illustration purposes, a simple 
implicit finite difference energy equation which includes 
a heat source term, Q v is given by equation (5). The 
heat source term, Qj, is treated as a latent heat term 
during the m At m+1 time step for the ith node. The 
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Figure 12. — Solidification initiation logic. 


superscript k’s represent the iteration during computation 
of latent heat contribution over the m At m+1 time step. 
An initial guess (k=0) of T i m+1 is computed from the 
solution routine such as equation (5). During the initial 
(k=0) estimate of T- m+1 , the source term, °Q. is equal 
to zero. A heat source term for the next (k+1) iteration, 
k+1 Q i? is computed from equation (6) by using the k Tj m+1 
value calculated from equation (5). If the convergence 
criterion, equation (7), is not satisfied, the k+1 Qi is 
used in the subsequent iteration to compute the new 
(k+1) estimate of T im+1 . The convergence of Q i 
correctly requires Tj m+1 — > T lat , which is observed in 
equation (6). After equation (7) has been satisfied, the 
initial latent heat contribution, Q i lat , is set equal to 
The initial liquid mass fraction at t=t m+1 is then com- 
puted by using equation (2) or (3). 



_ C,(T„ - %,,) ^ 

(6) 

At 

' Vi 

k+1 Qi,os - k Qi,os 

< QERR 

(7) 
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Figure 13— Melting completion logic. 


The iterative scheme to initiate melting of a given 
ith node was implemented in the thermal model by 
using the overall approach diagrammed in figure 11. 
Similarly, the logic that was developed to initiate 
solidification is shown in figure 12. The actual vari- 
able names used in the model are also used in the 
flowcharts (figs. 9 to 14). The flowchart variable names 
T, Told, Q, and Qold in figures 11 and 12 correspond 
t0 T i,m+l> T i,m> k+I Qi> andk Qi- When the flag NPHASE 
is equal to -3 or 3, initiation of melting or freezing 
occurs, and the iterative scheme described above is 


executed. Iteration limits are also observed in fig- 
ures 11 and 12. If the convergence criteria are not 
satisfied after 500 iterations, the simulation is stopped 
and a diagnostic message is written to the log file. 

Criteria must be provided to determine when phase 
change of the ith node is completed. Similar to the 
phase change initiation problem, the heat transferred to 
the node results in both latent and sensible heating 
during the n At n+1 time step. Analogous to the initia- 
tion of phase change case, an iterative scheme is again 
developed. The objective, however, is to determine the 
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Figure 14. — Solidification completion logic. 


sensible heat contribution and resulting temperature 
rise in the n At n+1 time step. 

The completion of phase change time step is defined 
as n At n+1 = t n+1 - t n . The times t n and t n+1 are 
determined by satisfying the criteria given in equa- 
tion (8) or (9) for melting or solidification, respectively. 

X i>n < 1 and X i>n+1 > 1 (Melting) (8) 

X i n > 0 and X i>n+1 < 0 (Solidification ) (9) 


The general approach that was used to determine the 
latent heat contribution during initiation of phase change 
is used to compute the sensible heat contribution dur- 
ing completion of phase change. The excess latent 
heat, Q exs at the n At n+1 time step is computed from 
equation (10) or (11). On the initial pass (k=0 itera- 
tion) through the phase change completion logic, Q exs j 
is used to make two guesses of Tj n+1 . One estimate 
of T in+1 is from equation (12). The second estimate 
is obtained from the solution routine of equation (5), by 
using Q exs as the heat source term The value of 
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Ti n+1 obtained from equation (12) is defined as °T i n+1 , 
while the T in+1 value computed from equation (5) is 
1 T n+1 . A change in the heat source term, AQ V which 
is computed from equation (13), is then used to deter- 
mine the heat source for the next iteration, k+1 Q^. The 
phase change completion logic used in the model is 
shown in figures 13 and 14. The variable names Hlat, 
m, TITR, Q', QCHNG, Q, and NITR used in both the 
flowcharts and the model correspond to h lat> m i> T est,i* 
Qexs,i> A Qi’ Qi> and k - 

hlat m i(l- Xi) 

Qexs.i = (Melting ) (10) 

At 

h lal mi(Xi) 

Q exs i = — - (Solidification) (11) 

’ At 

Testa = + Tlat (12) 


AQi = 


C l k + lnp k>r> \ 

i ( M,n+1 “ M,n+1 j 

At 


(13) 


k+1 Qi ■ 

= k Qi - AQi 

(14) 

AQj 

< QERR 

(15) 


A QERR value of 1 J/kg was used for the thermal model 
simulations discussed later. The maximum nodal latent 
heat content error for the smallest PCM nodes is 0. 1 percent, 
while the maximum error in nodal latent heat content 
for the largest PCM node is 0.05 percent. The nodal 
latent heat content error, t|, is defined as 


QERR 

mi hlat 


(16) 


Simulation of Combined Heating and Cooling Cycles 

Slight differences in model geometry exist between 
the heating and cooling modes. The differences result 
from the open and closed positions of the MLI shutter 
during cooling and heating, respectively. The approach 
that is used to address the model geometry differences 
relies heavily on various features of the thermal ana- 
lyzer, SINDA/FLUINT, that was chosen for the analysis. 


2 Cullimore, B.A., et al.: SINDA '85/FLUINT User's 
Manual, Version 2.3. NASA Contract NAS9-17448; 
MCR-90-512, Martin Marietta Corporation, Denver, 
CO, 1990. 


Therefore, the terminology and nomenclature used in 
this section requires some familiarity with SINDA/ 
FLUINT. 2 

During the melt/freeze cycle simulations, the open 
and closed positions of the radiation shutter must be 
modeled to properly simulate cooling and heating por- 
tions of the cycles, respectively. The radiator flare 
radiates directly to the GAS lid during the cooling 
portion of a cycle, while the radiation shutter is simu- 
lated closed for the heating portion of the cycle. The 
arithmetic node 116 (see figure 5) represents the radia- 
tion shutter in the closed position. Two sets of radia- 
tion conductances, one for the heating mode (Heating 
RADK’s), and another set for the cooling mode (Cool- 
ing RADK’s) are generated for the model. For the 
closed shutter (heating) mode, the Heating RADK’s are 
computed to the outside MLI surfaces (from the MLI to 
the GAS) and to the inside surfaces (from the TES 
specimen to the MLI) using two separate TRASYS models. 
The Cooling RADK’s consist of radiation conductors 
from the TES specimen to the surrounding MLI, radia- 
tion flare to the GAS lid, and radiation conductors from 
the MLI to the GAS. 

Both sets of radiation conductances cannot be applied 
simultaneously. The Cooling RADK’s must be “turned 
off’ during the heating mode, while the Heating RADK’s 
must be “turned off” during the cooling mode. There- 
fore, the actual conductance values of the Cooling RADK’s 
are saved in a dummy array during the heating mode. 
The Cooling RADK’s in the conductance block of the 
model are then set to zero. During cooling, the actual 
values of the Cooling RADK’s are restored. The con- 
ductance values of the Heating RADK’s are then trans- 
ferred to a dummy array and zeroed out during the 
cooling mode. 

The transfer of conductance values to and from the 
dummy arrays occurs in the operations block. After 
completion of the melting or freezing criteria in vari- 
ables 2, TIMEND is set equal to TIMEO. This stops 
execution of the network solution subroutine (in this 
case FWDBCK) and returns control to the operations 
block. By use of a DO LOOP which encloses the 
network solution subroutine, the melt/freeze cycle sce- 
nario can be simulated for a user specified number of 
orbits. 


Thermal Model Simulations 

The data presented and discussed herein are intended 
principally for comparison with other computer pro- 
grams that have recently been developed to predict the 
behavior of PCM with large volume changes (refs. 2 to 
4). Combined conduction and constant density phase 
change capabilities of the recetnly developed computer 
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programs can be benchmarked with these results. This 
section also serves to document the initial thermal modeling 
that was performed in support of the TES-1 experiment 
design effort. The modeling approach was strongly 
influenced by the need to provide conservative thermal 
analysis results for design purposes. Predictions of 
power requirements, peak temperatures, and freezing/ 
melting times were performed as part of the TES-1 
conceptual design. Model approximations, assumptions, 
and limitations that were influenced by design consid- 
erations are briefly summarized below. 

As indicated in the section titled Phase Change Pro- 
cedure and Logic, the SINDA thermal model is not 
capable of predicting the actual phenomena occurring 
in the TES material. Phenomena which cannot be simulated 
directly with the thermal model include the translucent 
behavior of the lithium fluoride, the 30 percent volume 
change of the salt, and the buoyancy driven and Marangoni 
flows. Both the lack of required property data and the 
complexity of the semitransparent behavior require major 
simplifying assumptions for the present thermal analy- 
sis. Void formation also changes the thermal characteristics 
of the TES canister, because the proportion of heat 
transfer from the canister wall to the TES material by 
radiation and conduction changes as the void grows 
and shrinks. Therefore, conservative assumptions were 
applied, when possible, that would predict greater total 
power and higher canister wall temperatures than would 
actually occur during the experiment. 

The following approximations and assumptions should 
be kept in mind while reviewing the results or compar- 
ing them with other data: 

(1) Buoyancy driven and Marangoni flows are neglected. 
Therefore, the thermal resistance of the TES material is 
greater, and this requires more heating and cooling 
energy to melt and solidify the TES material, respec- 
tively. The canister wall peak temperatures also increase 
(conservative). 

(2) The void is placed in the realistically worse location 
for heating. The void is expected to migrate to a hot 
spot. For the TES-1 experiment, this will be at the 
outer wall. Placing the void along the whole outer 
canister wall requires that heat be transferred to the 
TES material by radiation. A small amount of heat is 
also conducted along the wall and then is conducted 
from the top and bottom canister walls into the TES 
material. This assumption increases both the peak can- 
ister wall temperature and the energy required to melt 
the TES material. 

(3) Heat losses from thermocouples are not included 
in the simulation. Conduction losses from the heater 
through the heater base are also neglected 
(nonconservative). 

The following is a description of the characteristics 
and scenarios modeled in each of the four simulations 


identified in table V. Each simulation is divided into 
two scenarios-an initial heatup of the TES specimen, 
and the TES melt/freeze cycles. The initial heatup 
simulation begins with the initial temperatures of 300 K 
and continues through complete melting of the PCM. 
This simulates the start of the actual experiment on 
board the orbiter. The effects of a prolonged cooldown 
period prior to initiation of the experiment are neglected. 
However, these cooldown periods for different orbiter 
orientations (e.g., Sun facing and Earth viewing) could 
be incorporated in subsequent simulations of the flight 
experiment. 


TABLE V.-TES EXPERIMENT 
SIMULATIONS 


Simulation 

run 

Lid mass, 
kg 

PCM, K|jq, 
W/m-K 

1 

20.34 

1.73 

2 

9.09 

1.73 

3 

34.09 

1.73 

4 

20.34 

3.46 


The first simulation in table V (run 1), is defined as 
the base case. The TES-1 experiment is placed in a 
conventional GAS can with no modifications. Esti- 
mates of power, peak temperatures, and freezing/melt- 
ing times were established from run 1. Runs 2 and 3 
differ from run 1 only in the GAS lid mass. Results 
from runs 1, 2, and 3 were used to determine the impact 
of the GAS lid mass on melting and solidification times. 
The fourth simulation (run 4) uses a LiF liquid thermal 
conductivity that is twice the value of the base case 
(run 1). According to reference 8, combined radiation 
and conduction through the liquid LiF can be approxi- 
mated by using a thermal conductivity twice the value 
of K liq . Run 4 was conducted to investigate the TES-1 
experiment performance with the use of a more realis- 
tic approximation of the heat transfer through the PCM 
than the conservative approximation used for run 1. 
The solid phase of LiF is also translucent; however, 
this was not included in the run 4. 

Both the initial heatup and the melt/freeze cycle sce- 
narios remain the same for all four simulations. During 
the initial heat-up, the radiation shutter above the TES 
radiation flare remains closed. The heater is energized 
with a constant power of 260 W. The environmental 
temperature which is seen by the external surfaces of 
the GAS and the radiation shutter is 250 K. Typical 
environmental (sink) temperatures for an Earth view 
orbit are on the order of 250, while sink temperatures 
to deep space are generally less than 100 K. Results of 
initial heat-up simulation (heat-up time, heat-up energy, 
and TES specimen temperatures) are insensitive for 
environmental temperatures less than 300 K. Results 
of primary interest from the initial heatup scenario are 
the time to begin melting, the time of complete melting, 
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peak TES canister wall temperatures, heater power, and 
required heater energy. 

The final temperatures from the initial heatup sce- 
nario were used as initial temperatures for the melt/ 
freeze cycle simulation. Four and one half cycles, 
beginning with the solidification of the LiF and ending 
with the solidification of the LiF, were simulated. After 
complete solidification of the lithium fluoride, the melting 
portion of a given cycle was initiated. The melting 
process was completed after complete melting of the 
LiF. Details of the combined heating and cooling cycle 
simulation were presented earlier. The environmental, 
or sink, temperature during the melt/freeze simulation 
was 250 K. The external surface of the GAS lid and all 
external sides of the GAS can radiate to space (envi- 
ronment). The TES radiator flare radiates to the GAS 
lid during the PCM solidification. Simulation results 
are again not affected by sink temperatures of 300 K or 
less. All heat losses from the GAS are assumed to be 
by radiation to the environmental (sink) temperature. 


MLI Q 
MLI Pavg 


Base Q 


Base Pavg 


U TES 


U Dot 


Q Lat 


MLI heat loss during heating 

average heat-loss power through MLI 
during heating portion of given cycle 
(MLI Q/Delta Time/60) 

heat loss from TES mounting base to 
node 109 

average heat-loss power through mount- 
ing base for corresponding Delta Time 

PCM internal energy change for cor- 
responding Delta Time 

average cooling or heating rate of PCM 
for corresponding Delta Time 

PCM latent heat change for corre- 
sponding Delta Time 


Q Sense 

Results and Discussion 

Results consist of tabulated results provided in tables VI 
to IX, transient temperature profiles in appendix B, and T Lid 
color temperature contour plots in appendix C. The 
following nomenclature and definitions describe the 
results presented in tables VI to IX: 

T Htr 

Compl Time completion time from beginning of 
scenario 


TES specimen (nodes 1 to 58) sen- 
sible heat change for corresponding 
Delta Time 

temperature of the GAS lid (node 108) 
at corresponding time into simula- 
tion, Compl Time 

peak temperature of heater nodes at 
corresponding time into simulation, 
Compl Time 


Delta Time melting/freezing time per cycle 
Heater Q heater energy for given cycle 


T Canister -OD peak temperature of outer TES canis- 
ter wall at corresponding into simu- 
lation, Compl Time 


TABLE VI.— RUN 1 PHASE CHANGE COMPLETION RESULTS 
[Lid mass, 20.3 kg.] 



Initial 

heatup 

Cycle 1 

Cycle ; 

l 

Cycle 

3 

Cycle 4 

Cycle 5 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Compl time, sec 

7993.0 

2015.0 

5864.0 

7872.0 

11774.0 

13773.0 

17674.0 

19674.0 

23573.0 

25573.0 

Delta time, min 

133.0 

33.6 

64.2 

33.5 

65.0 

33.3 

65.0 

33.3 

65.0 

33.3 

Hpotpr O kJ 

2078.3 


1001.0 


1015.0 


1014.0 


i n i a n 


I H'dlCI y , AJ 

MI I O kJ 

-236.0 


-202.2 


-218.0 


-218.0 


1 U 14. u 
o i o n 


MI I Pavp W 

-29.6 


-52.5 


-55.8 


-55.8 


-Z 1 o .u 

-55.8 


IVlJLtl I a v TT 






Base Q, kJ 

-133.0 

-42.3 

-67.1 

-42.1 

-61 A 

-41.8 

-67.3 

-41.8 

-67.3 

-41.8 

Base Pavg, W 

-16.6 

-21.0 

-17.4 

-21.0 

-17.3 

-20.9 

-17.3 

-20.9 

-17.3 

-20.9 

U TES, J 

1709.0 

-736.0 

731.4 

-731.4 

729.2 

-728.6 

729.0 

-729.0 

728.9 

-729.0 

U Dot, W 

213.8 

-365.0 

190.0 

-364.0 

186.9 

-364.5 

186.9 

-364.4 

186.9 

-364.4 

Q Lat, kJ 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

Q Sense, kJ 

1388.0 

-414.6 

409.9 

-409.9 

407.7 

-407.2 

407.6 

-407.3 

407.4 

-407.4 

T Lid, K 

281.0 

295.0 

285.0 

298.0 

286.2 

299.0 

287.0 

299.0 

288.0 

300.0 

T Htr, K 

1276.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

T Canister-OD K 

1237.0 

1065.0 

1234.0 

1065.0 

1234.0 

1065.0 

1234.0 

1065.0 

1234.0 

1065.0 
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TABLE VII.— RUN 2 PHASE CHANGE COMPLETION RESULTS 


[Lid mass, 9.1 kg.] 



Initial 

Cycle 1 

Cycle : 

2 

Cycle 

3 

Cycle 4 

Cycle 5 


heatup 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Compl time, sec 

8045.0 

2005.0 

5818.0 

7826.0 

11727.0 

13727.0 

17627.0 

19628.0 

23528.0 

25528.0 

Delta time, min 
Heater Q, kJ 

134.0 

2091.7 

33.4 

63.6 

991.4 

33.5 

65.0 

1014.3 

33.3 

65.0 

1014.0 

33.4 

65.0 

1014.0 

33.3 






MLI Q, kJ 

-252.0 

-194.2 

-217.7 

-217.7 

-217.8 






MLI Pavg, W 

-31.4 

-50.9 

-55.8 

-55.8 

-55.9 






Base Q, kJ 

-133.6 

-42.1 

-66.5 

-42.1 

-67.4 

-41.8 

-67.3 

-41.8 

-67.3 

-41.8 

Base Pavg, W 

-16.6 

-21.0 

-17.4 

-21.0 

-17.3 

-20.9 

-17.3 

-20.9 

-17.3 

-20.9 

U TES, J 

1705.6 

-732.0 

730.8 

-731.0 

729.1 

-728.7 

729.0 

-728.9 

728.9 

-728.6 

U Dot, W 

212.0 

-36.5 

191.7 

-36.4 

186.9 

-364.3 

186.9 

-364.3 

186.9 

-364.3 

Q Lat, kJ 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

Q Sense, kJ 

1384.1 

-410.6 

409.3 

-409.5 

407.7 

-407.2 

407.5 

-407.5 

407.4 

-407.1 

T Lid, K 

276.0 

298.1 

282.0 

301.5 

283.1 

302.7 

283.9 

303.2 

284.3 

303.4 

T Htr, K 

1275.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

T Canister-OD K 

1236.0 

1065.0 

1234.0 

1065.0 

1234.0 

1065.0 

1234.0 

1065.0 

1234.0 

1065.0 


TABLE VIII. — RUN 3 PHASE CHANGE COMPLETION RESULTS 


[Lid mass, 34.1 kg.] 



Initial 

heatup 

Cycle 1 

Cycle : 

2 

Cycle 

3 

Cycle 4 

Cycle 5 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Compl time, sec 

8043.0 

2005.0 

5906.0 

7905.0 

11806.0 

13806.0 

17707.0 

19706.0 

23606.0 

25606.0 

Delta time, min 

134.0 

33.4 

65.0 

33.3 

65.0 

33.3 

65.0 

33.3 

65.0 

33.3 

Heater O kJ 

2091.2 


1014.0 


1014.3 


1014.3 


inii n 


Ml I O kJ 

-252.3 


-218.0 


-217.7 


-217.9 


1U14.U 


i VI JLs 1 y , IVJ 

MI I Pave W 

-31.4 


-55.8 


-55.8 


-55.8 


-Z 1 / .0 

-55.8 


1V1JL1 r CL ▼ ▼▼ 






Base Q, kJ 

-133.5 

-42.1 

-67.5 

-41.8 

-67.4 

-41.8 

-67.3 

-41.8 

-67.3 

-41.8 

Base Pavg, W 

-16.6 

-21.0 

-17.3 

-21.0 

-17.3 

-21.0 

-17.3 

-20.9 

-17.3 

-20.9 

U TES, J 

1705.5 

-7320.4 

729.1 

-728.8 

729.1 

-729.1 

729.1 

-728.7 

729.1 

-729.0 

U Dot, W 

212.0 

-36.5 

186.9 

-364.6 

186.9 

-364.6 

186.9 

-364.5 

186.9 

-364.5 

Q Lat, kJ 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

Q Sense, kJ 

1384.0 

-410.6 

407.6 

-407.4 

407.7 

-407.6 

407.6 

-407.3 

407.6 

-407.5 

T Lid, K 

284.8 

294.0 

286.8 

295.7 

288.0 

296.7 

288.8 

297.4 

289.3 

297.8 

T Htr, K 

1275.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

1274.0 

1026.0 

T Canister-OD K 

1236.0 

1065.0 

1234.0 

1065.1 

1234.0 

1065.0 

1234.0 

1065.0 

1234.0 

1065.0 


All results in column 1 of the tables represent the final 
values from the initial heat-up scenario. All subse- 
quent columns contain results from the melt/freeze cycle 
scenarios. 

A comprehensive set of transient temperature pro- 
files from run 1 is given in appendix B. Run 1 tempera- 
ture profiles are shown for all the LiF PCM nodes 
(figs. B1 to B7); the TES canister inner and outer wall 
nodes 64 to 73 (figs. B8 and B9); heater nodes 97, 99, 
101, and 103 (fig. BIO); radiator surface nodes 11, 12, 
14, 16, and 18 (fig. B 1 1 ); conductor-rod nodes 11, 54, 
56, and 58 (fig. B12); lid and GAS wall nodes 105 to 
108 (fig. B13); mounting-base nodes 110, 111, and 1 12 
(fig. B14); and spacer nodes 109 and 9109 (fig. B15). 
Figures B16 and B17 provide a comparison of temperature 


profiles for runs 1, 2, and 3 of LiF PCM node 46, and 
the GAS lid node 108, respectively. 

Six temperature contour plots are given in appendix C. 
Temperature contours of the TES specimen after completion 
of the initial heatup scenario and after the melt/freeze 
cycles scenario are shown in figures C.l and C.2, 
respectively. Temperature contours for the overall GAS 
configuration are similarly shown in figures C.3 and 
C.4. Figure C.5 shows eight temperature contours 
from run 1 during the cycle-3 freezing and thawing 
process. The contours represent (a) the approximate 
first quarter of the freeze cycle, (b) the approximate 
midpoint of the freeze cycle, (c) the approximate third 
quarter of the freeze cycle, (d) complete solidification, 
(e) the approximate first quarter of the melt cycle, (f) the 
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TABLE IX. — RUN 4 PHASE CHANGE COMPLETION RESULTS 


[Thermal conductivity, 2x Liquid Kjjp.] 



Initial 

heatup 

Cycle 1 

Cycle 2 

Cycle 3 

Cycle 4 

Cycle 5 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Melting 

Solidification 

Compl time, sec 

7878.0 

1922.0 

5670.0 

7597.0 

11343.0 

13269.0 

17016.0 

18943.0 

22688.0 

24615.0 

Delta time, min 

131.3 

32.0 

62.5 

32.1 

62.4 

32.1 

62.5 

132.1 

62.4 

32.1 

Heater Q, kJ 

2048.3 


974.5 


974.0 


974.0 


Q 70 7 


MLI O kJ 

-236.5 


-202.8 


-202.6 


-202.9 


y / j. / 


MLI Pavg, W 

-30.0 


-54.1 


-54.1 


-54.2 


-ZUZ . J 

-54.1 








Base Q, kJ 

-131.6 

-40.5 

-65.9 

-40.5 

-65.7 

-40.4 

-65.7 

-40.4 

-65.6 

-40.4 

Base Pavg, W 

-16.7 

-21.1 

-17.6 

-21.0 

-17.5 

-21.0 

-17.5 

-21.0 

-17.5 

-21.0 

U TES, J 

1680.3 

-705.7 

705.8 

-705.8 

705.7 

-705.4 

705.6 

-705.6 

705.5 

-705.5 

U Dot, W 

21.2 

-367.0 

188.3 

-366.3 

188.4 

-366.2 

188.3 

-366.2 

188.4 

-366.1 

Q Lat, kJ 

213.3 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

321.5 

-321.5 

Q Sense, kJ 

1358.8 

-384.3 

384.3 

-384.3 

384.2 

-383.9 

384.2 

-384.1 

384.1 

-384.1 

T Lid, K 

281.0 

294.4 

284.5 

297.1 
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approximate midpoint of the melt cycle, (g) the approximate 
third quarter of the melt cycle, and (h) complete melt. 
A temperature contour of the initial heatup scenario 
from run 4 is given in figure C.6 for comparison with 
the base case in figure C.l. 

Base-case values of melting time, peak temperatures, 
and PCM heat changes for the initial heatup scenario 
are given in column 1 of table VI. During the initial 
heatup simulation, melting begins after approximately 
1.5 hr. Complete melting occurs after 2.17 hr from the 
beginning of the initial heatup time. The peak canister 
wall temperature of 1237 K is less than the 1300 K 
limit imposed during the conceptual design of the TES 
experiment. This temperature is also the maximum 
temperature that occurs over both scenarios, the initial 
heatup and melt/freeze cycles. The TES specimen tem- 
perature contours in figure C3 indicate that the domi- 
nant direction of heat transfer is the radial direction. 
However, the heat losses from both the radiator flare 
and the mounting base which are evident in figure C3 
also result in axial heat flow. Inspection of run 2 and 
run 3 initial heatup results in tables VII and VIII show 
negligible differences from the run 1 results. 

The base-case heat loss from the TES specimen by 
radiation through the MLI and by conduction through 
the base is 366 kJ (the sum of rows 4 and 6 in table VI) 
for the initial heatup scenario. The radiative losses 
from the TES specimen are the heat losses from the 
MLI surrounding the TES specimen and the losses through 
the closed radiation shutter above the TES specimen. 
The average rate of heat loss (power) from the TES 
specimen through complete melting is 46 W. During 
the initial heatup of the TES specimen, the GAS 
temperatures decrease as shown in figure B13. The 
GAS wall temperatures are approximately 280 K at the 
completion of the initial heatup scenario. The tempera- 


ture contour plot, figure C3, shows uniform GAS wall 
temperatures at the completion of the initial heatup 
while, temperature gradients exist throughout the TES 
specimen and mounting assembly. 

Nominal cycle values of melting/freezing cycle times, 
temperatures, and PCM heat changes are given in table X. 
These values are obtained from the run 1 cycle results; 
however, they are also typical of runs 2 and 3. Melt 
and freeze times of 65 and 33 min are within acceptable 
limits of 1 hr +15 min -6 min and 34 min + 6 min -4 min 
as specified in the TES Experiment Technical Require- 
ments Document (TRD) (A NASA Lewis Internal Docu- 
ment prepared by David Namkoong in 1990). 

TABLE X. -PREDICTED NOMINAL TES 
PERFORMANCE RESULTS PER THE 
MELT/FREEZE CYCLE 


Melting time, min 65 

Freezing time, min 34 

Peak canister temperature, K 1234 

Peak heater temperature, K 1274 

Heater energy, kJ 1014 

PCM latent heat change, kJ 322 

TES specimen sensible heat change, kJ 407 


During the melt/freeze cycles, the base-case canister- 
wall peak temperature is 1234 K, which is less than the 
peak temperature at the completion of the initial heatup 
scenario. Transient temperature profiles of LiF nodes 
are given in appendix B, figures B1 and B2. Specific 
node locations can be found from the node diagrams 
(figs. 5 and 7). The LiF PCM transient profiles in 
figures B1 to B5 show maximum subcooling of the 
PCM of approximately 90 K below the melting tem- 
perature of 1121.3 K. During heating, a temperature 
rise of approximately 85 K above melt temperature is 
observed for the PCM. The temperature contour plot 
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for the completion of cycle-5 solidification shows coolest 
temperatures at the radiator flare surface as expected. 
In addition, the temperature profile of the TES speci- 
men is observed to be strongly two-dimensional. Mounting- 
base temperature gradients seen in the overall GAS 
configuration contours (figs. C3 and C4) result in the 
base conduction losses shown in tables VI to IX. 

The base-case temperature contour plots in figure C5, 
show the freeze/melt process through the complete cycle 3. 
As solidification occurs, the LiF is observed to freeze 
inward from all TES canister walls. The lower central 
canister region is the last to freeze. It may be desirable 
to require LiF nearest the outer wall to freeze last, as 
suggested in the 1990 NASA Lewis Internal TES Tech- 
nical Requirements Document by David Nankoong. 
Therefore, applying a small amount of power to the 
heater may be necessary during the solidification pro- 
cess. During the early stages of melting, as shown in 
figures C5(e) and (f), heat flow in the PCM is primarily 
in the radial direction. At the completion of melting 
(fig. C5(h)) the lower portions of the conductor rod and 
PCM are observed to have cooler temperatures. These 
cooler temperatures result from the conduction losses 
to the mounting base. 

The TES specimen heat losses through the MLI and 
through the base during the heating portion of a cycle 
were 218 and 67 kJ, respectively. The resulting aver- 
age heat loss rates were 56 and 17 W for the MLI and 
base, respectively. Recall that an estimated or desired 
value of the average MLI heat loss rate is used to 
determine the MLI effectiveness as discussed in the 
Model Geometry and Thermal Properties section. Summing 
the TES internal energy change and heat losses (rows 
4, 6, and 8 in table VI) during the heating portion of the 
third cycle yields 1003 kJ. The resulting variation with 
the heater energy input is small (1.1 percent). 

A comparison of tables VI to VIII reveals that the 
only noticeable difference among runs 1 , 2, and 3 is the 
lid temperature during the melt/freeze scenario. The 
PCM temperature profiles for LiF node 46 are plotted 
in figure B16 for runs 1, 2, and 3. The maximum and 
minimum temperatures of the three runs are the same. 
However, the melt/freeze cycle times of the smallest 
mass GAS lid (run 2) lag the simulations with larger 
GAS lid masses. The lagging of the run-2 temperature 
profile increases in each cycle. Figure B17 shows a 
comparison of transient temperature profiles for the 
GAS lid, node 108, from runs 1, 2, and 3. Run 2, which 
has the smallest GAS lid mass (9.1 kg), shows the larg- 
est temperature swings and the greatest temperature 
rise with each cycle. The increasing lag of the run-2 
temperature profile in figure B16 may be caused by the 
gradual increase in GAS lid temperature in subsequent 
cycles. An increase in lid temperature also causes an 
increase in freezing time. Run 3, which has the largest 


TABLE XI. — GAS LID TEMPERATURE RISE 
AS FUNCTION OF LID MASS 


Simulation 

run 

Lid mass, 
kg (lb) 

Lid temperature 
rise, 

K 

1 

20.34 (44.25) 

19.0 

2 

9.091 (20) 

27.4 

3 

34.09 (75) 

13.0 


GAS lid mass (34.1 kg), shows the smallest tempera- 
ture swing and smallest temperature rise per cycle. 
Table XI gives the temperature rise of the GAS lid 
during the 4.5 freeze/melt cycles for different GAS lid 
masses (simulations). 

Run 4 represents a more realistic simulation of the 
TES experiment by attempting to account for heat trans- 
ferred to the LiF by radiation in addition to the heat 
transferred by conduction. An approximation to account 
for the radiative contribution is suggested in reference 8. 
This approximation consists of simply doubling the 
liquid thermal conductivity value of the LiF. Compari- 
son of the run-4 simulation results with those of the 
base case, run 1, reveals that the melting and freezing 
times will decrease by approximately 2.5 and 1.2 min, 
respectively. The canister-wall peak temperature of 
1215 is 22 K lower than that of the base case. These 
comparisons translate into a 3.2-percent decrease in the 
TES specimen internal energy change and a 4-percent 
decrease in heater energy per cycle with respect to the 
base case. Figure C6 shows the cooler temperatures of 
the heater, canister walls, conductor rod, and radiator 
flare at the completion of the initial heatup from run 4. 


Concluding Remarks 

An approach for modeling phase change using ther- 
mal analyzers was successfully developed and imple- 
mented. The advantages of this approach with respect 
to the effective Cp approach are insensitivity to computer 
roundoff error and ability to maintain a constant phase 
change temperature. The method is applicable to ther- 
mal analyzers that solve the energy equation in terms 
of temperature (this includes most commercial analyz- 
ers). The method can be incorporated into thermal 
analyzers, such as SINDA, which permit access to basic 
solution variables (i.e., temperatures, thermal capaci- 
tances, and thermal conductances). 

The thermal model was developed to provide conser- 
vative estimates of required heater energy, melting and 
freezing times, and peak temperatures. Conservative, 
for this analysis, implies greater predicted heater energy, 
larger predicted times and hotter predicted peak tem- 
peratures with respect to actual values. For all simu- 
lations, the cycles were repeatable after the first melt/ 
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freeze cycle, and the PCM freezing and melting times 
were within acceptable limits. Peak temperatures of 
the PCM canister did not exceed the 1300 K tempera- 
ture limit. The aluminum GAS wall temperatures never 
exceeded 300 K, which is well below GAS or elec- 
tronic components overheat (failure) temperature lim- 
its. Three simulations were performed to investigate 
the effects of the GAS lid mass on the melting and 
freezing times of the PCM. Results from the lid mass 
comparison indicate negligible effect on the freezing 
and melting times for GAS lid masses from 9.1 to 
34.1 kg. Therefore, the use of a lid mass between 9.1 
and 34.1 kg will not affect the TES experimental results. 
The fourth simulation uses a simple approximation to 
account for combined radiation and conduction trans- 
port through the liquid LiF (the semitransparent nature 
of the solid phase was not included). Comparison of 
the simulation results with those of the base case, run 1, 
reveals that the melting and freezing times will decrease 
by approximately 2.5 and 1.2 min, respectively, while 
the peak canister wall temperature is 22 K lower than 
the base case. The liquid semitransparent approximation 
results in 3.9-, 3.8-, and 2-percent decreases in freezing 
times, melting times, and peak canister temperatures 
with respect to conduction only. The TES specimen 
internal energy change and heater energy per cycle 
decrease by 3.2 and 4 percent, respectively, compared 
to the base case. For all simulations, the influence of 
buoyance or surface tension driven convection is expected 
to reduce melting/freezing times and peak PCM canis- 
ter temperatures. 

The thermal model and results can be used to evalu- 
ate (benchmark) new computer programs or models that 


have been developed to predict solid/liquid phase change 
behavior. Specifically, results are provided for 
benchmarking combined conduction and phase change 
capabilities of the unverified computer programs or 
models. 
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Appendix A 


Finite Difference Solutions and 
Development of Phase Change Initiation Logic 


Three finite difference forms of the energy equa- 
tion - explicit, implicit, and explicit/implicit (Crank- 
Nicholson approximation) - are given in equations (Al) 
to (A3). All three equations consist of the sensible heat 
storage term, a heat-source term, and conduction terms. 
For the present discussion, other possible terms such as 
radiation conductance terms are neglected. 


Ci(T iin +i - T i>n ) 
At 


= Qi 


+ 




(explicit) 


(Al) 


t n+1 , can be computed immediately and independently. 
From inspection of equations (A2) and (A3), we note 
that Tj n+1 is also dependent on the other node tempera- 
tures, Tj’s at time, t^. The implicit and Crank-Nicholson 
solution schemes are attractive because they are uncon- 
ditionally stable (not necessarily correct). 

The CG min , Cj/XGjp is an indicator of the maximum 
stable time step size for an explicit solution for the TES 
thermal model the CG min varies approximately two 
orders of magnitude. Therefore, the Crank-Nicholson 
solution routine (SINDA/FLUINT Subroutine FWDBCK) 
was chosen for this analysis. The exact finite differ- 
ence equation that is solved in subroutine FWDBCK is 
given in equation (A5). 3 


^i/i) 

At 


= Qi 


+ 



- T u + i)] 


(implicit) 


(A 2) 



At 


= Qi 





(Crank - Nicholson ) 


(A3) 


Rearranging equation (Al) to solve for Tj n+1 leads to 
equation (A4). 



From equation (A4) we see that the temperature of an 
ith node at the new time value, t n+1 , can be “explicitly” 
solved for in terms of node temperatures at the old 
time, t n . Therefore, any temperature at the new time, 


3 Cullimore, B.A., et al.: S1NDA '85/FLUINT User's 
Manual, Version 2.3. NASA Contract NAS9- 17448; 
MCR-90-512, Martin Marietta Corporation, Denver, 
CO, 1990. 



The latent heat contribution, Q i? in the m At m+1 time 
step is exactly equal to Qj os (computed from eq. (4)) 
for an explicit solution routine. The equivalence of Qj 
and Qj os for an explicit finite difference solution is 
demonstrated below. 

First, we rewrite the explicit finite difference 
equation (Al) in a form given by equation (A6), which 
will allow us to apply the phase-change initiation scheme 
described in the section titled Simulation of Combined 
Heating and Cooling Cycles. The superscript k denotes 
the kth iteration at the m At m+1 time step. 





At 



(A 6) 


The remaining term on the right side of equation(A6) is 
invariant with k. Therefore, we can equate the left side 
of equation (A6) for the 0th and 1st iterations as shown 
in equation (A7). 
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Q(°Tj,n-H ' Tjji) 0 _ ^i( T i>n-H " ^ 


°Qi 


'*») 


Qi 


(A 7) 


At ' A At 

The temperature °T i n+1 by definition is also the over- 
shoot temperature, Tj os , while °Q l = 0. Substituting 
°T- +1 an ^ °Qi into equation (6) gives equation (A8) 


for l Q { . 


l Q[ = 


Ci(Tiat - \os) 


At 


(A8) 


Substituting equation (A8) into equation (A7) and sim- 
plifying yields equation (A9). 

Ti,os - T i>n = (\ n+l - Tj n ) - (T lat - T i>os ) (A9) 

Both T ios and Tj n cancel out of equation (A9) 
leaving tfie result: 


l. 


A i,n+1 


= Ti 


lat 


Since 1 T i>n+1 = T lat , 2 Q i = which satisfies the 
convergence criteria given by equation (7). In summary, 


equation (4) immediately (no further iterations) gives 
the latent heat contribution in the m At m+1 time step for 
an explicit finite difference solution. 

By applying a similar argument to an implicit finite 
difference solution, equation (A2), is rewritten as 
equation (A 10). The temperatures on the right side of 
equation (A10) are computed at t=t n+1 . Therefore, 
superscript k is now required on these terms. 


CijX+i - T u ) 
At 



(A10) 


In contrast to equation ( A6), the right side of equation (A 1 0) 
not invariant with k. Since all temperatures at the 
time = t n+1 must be solved for simultaneously, any 
change in Tj n+1 affects all other temperatures, Tj n+1 ’s. 
Hence, an iterative scheme such as the one developed 
in the section on Simulation of Combined Heating and 
Cooling Cycles is required to determine the latent heat 
contribution during the m At m+1 time step (initiation of 
phase change). 
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Appendix B 


Transient Temperature Profiles for Selected TES Components and Nodes 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B1. — LiF nodes, row 1. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B2 LiF nodes, row 2. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B3.— LiF nodes, row 3. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B4. — LiF nodes, row 4. 


27 


Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B6. — LiF nodes, column 2. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B 7. — LiF nodes, column 3. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B8. — Canister inner-wall nodes. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B9. — Canister outer-wall nodes. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure BIO. — Heater nodes. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B1 1 . — Radiator surface nodes. 
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Nodal temperature, K Nodal temperature, K 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 

Figure B12. — Conductor-rod nodes at inner radius. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B13. — Lid and GAS wall nodes. 
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Nodal temperature, K Nodal temperature, K 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B14. — Mounting-base nodes. 


37 


Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 
Figure B15. — Spacer nodes. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 

Figure B16. — Effect of lid mass on PCM temperature for PCM node 46. 
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Nodal temperature, K Nodal temperature, 



(a) Initial heatup scenario. 



(b) Freeze/melt cycle scenario. 

Figure B17. — Effect of lid mass on lid temperature for lid node 108. 
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Appendix C 


Temperature Contour Plots 
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Temperature, 

K 
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1176 
1166 
1157 
1148 
1138 I 


Figure Cl . — Run 1 temperature contours of TES specimen after 
initial heatup. Heater power, 260 W; environmental (sink) tem- 
perature, 250 K. 



Figure C3. — Run 1 temperature contours of overall 
TES/GAS configuration after initial heatup. Heater 
power, 260 W; environmental (sink) temperature, 
250 K. 


Temperature, 

K 



Figure C2. — Run 1 temperature contours of TES specimen after 
cycle-5 solidification. Heater power, 260 W; environmental 
(sink) temperature, 250 K. 



Figure C4. — Run 1 temperature contours of overall 
TES/GAS configuration after cycle-5 solidification. 
Heater power, 260 W; environmental (sink) tem- 
perature, 250 K. 
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Temperature, 

K 



(d) Completion of freeze cycle. 
Time, 13773 sec. 


(e) Approximate first quarter of 
melt cycle. Time, 14773 sec. 


(f) Approximate midpoint of melt 
cycle. Time, 15773 sec. 



(g) Approximate third quarter of (h) Completion of melt cycle, 

melt cycle. Time, 16773 sec. Time, 17674 sec. 


Figure C5. — Run 1 temperature contours of TES specimen during cycle 3. 
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Temperature, 

K 
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Figure C6. — Run 4 temperature contours of TES specimen after 
initial heatup. Heater power, 260 W; environmental (sink) tem- 
perature, 250 K. 
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